Fluctuations in a diffusive medium with gain 
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We present a stochastic model for amplifying, diffusive media like, for instance, random lasers. 
Starting from a simple random-walk model, we derive a stochastic partial differential equation for 
the energy field with contains a multiplicative random-advection term yielding intermittency and 
power-law distributions of the field itself. Dimensional analysis indicate that such features are more 
likely to be observed for small enough samples and in lower spatial dimensions. 
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Wave transport in disordered media can be described 
as a multiple scattering process in which waves are ran- 
domly scattered by a large number of separate elements 
p|. To first approximation this gives rise to a diffusion 
process. A particularly interesting situation arises when 
gain is added to a random material. In optics this is 
realized, for instance, in the form of a suspension of 
micro particles with added laser dye or by grinding a 
laser crystal. If the total gain becomes larger then the 
losses, fluctuations grow and these systems exhibit a las- 
ing threshold J2[ .yielding the so called random laser (see 
e.g. Refs. [3, |4| and the bibliography therein). The 
complexity generated by the interplay between gain and 
disorder leads to intriguing connections with other fun- 
damental problems like Anderson localization [j| or the 
physics of glasses @. Besides their fundamental interest, 
random lasers are likely to have a technological impact 
as low-cost light sources. 

Similar situations occur in other branches of physics 
like neutron diffusion in fissile materials, or stochastic 
wave growth Q and acceleration of plasma particles Q 
As known, the competition between growth and propaga- 
tion or diffusion is also a basic mechanism of population 
dynamics and theoretical ecology @. 

Diffusive random lasing has been observed experimen- 
tally in various active random media, including pow- 



dicted [21] . Remarkably, those prediction were very re- 
cently confirmed experimentally (see also Ref. 23 1). 

In this work we present first a simple stochastic model 
that can be analytically solved and that yields large sta- 
tistical fluctuations. From it, we derive an equation for 
the energy field which contains a stochastic multiplica- 
tive term whose strength controls such anomalous fluc- 
tuations. The magnitude of this term introduces another 
scale in the problem which thus provides a criterion for 
the observability of such intermittency in the field distri- 
bution. 

We choose to describe isotropic diffusion of light in 
terms of an ensemble of independent random walkers 
each carrying a given energy (number of photons). This 
may be visualized as an ensemble of "beams" propagat- 
ing independently throughout the sample, each interact- 
ing with an underlying atomic population providing a 
gain mechanism via stimulated emission. This procedure 
of attaching an energy to a random walk of photons is a 
common practice in Monte Carlo simulations of absorp 
tion in complex materials like e.g 
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ders, laser dye suspensions and organic systems 10h12|. 
Theoretical descriptions often relies on the diffusive ap- 
proximation either by reaction-diffusion type of equa- 
tions @, [l3| or, on a more microscopic level, by the 



tissues 

been also employed previously for random lasers 
under the implicit assumption that all the photons gen- 
erated by amplification are diffused in the same direction 
at each scattering event. 

Let us denote by x and E the walker position and 
energy. We discuss the one-dimensional case in which 
the walker resides on a finite interval, < x < L. The 
dynamics is formulated as follows. A new walker is gen- 
erated at a random position by a spontaneous emission 
event with a rate 7 and initial energy E = e. In terms of 
the underlying active media, 7 denotes the spontaneous 
emission rate of the single atom. The walker position x 
changes to x ± a according to a standard random walk 
rule on a lattice with spacing a. At the same time, the 
walker energy E may increase by one unit due to the 
process of stimulated emission, E — > E + e, with a rate 
T(E). The simplest choice would be T(E) = or, to 
mimic saturation effects we may consider a gain of the 
form T(E) = jE/(l + E/E s ). 

The probability Pi >n for the walker to be at x = ia 
and crossovers among different statistics has been pre- having an energy E = ne evolves according to a master 



master-equation approach [14] . Monte Carlo simulations 
play of course an important role [ljl [l6| . 

One of the salient experimental features of random 
laser emission is its large statistical variability 



17]. In- 



deed, already within the diffusive approximation, the ad- 
dition of gain (and saturation) naturally generates fat- 
tailed distributions [18| that stems from rare long light 
paths [16] . This mechanism is of pure statistical origin 
and that does not require localization or interference [17] . 
It was also proposed that such systems can exhibit Levy 



type statistics in the distribution of intensities 18l420l| 
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equation, that can be solved by standard methods [25j . 
Here, to simplify further, we work directly in the con- 
tinuum limit and treat x and E as continuous variables, 
obeying the Langevin equations 



x = V2D£, E = T(E) + y/T(E)ri 



(1) 



where £, rj are ^-correlated Gaussian variables with unit 
variance. To keep things as simple as possible, bulk ab- 
sorption or the possibility that diffusion is affected by 
the energy are neglected from scratch but can be easily 
included. 

To demonstrate that this dynamics naturally gener- 
ate power-law distributions of energies, we solve the as- 
sociated Fokker-Planck equation (Ito interpretation of 
Eq. ©) 



P = D 



d 2 p 



c) 



dx 2 dE 



2 dE 



(2) 



where P(x, E, t) is the probability of finding a walker 
with energy E at x. Four boundary conditions arc nec- 
essary on the contour of the domain [0,L] x [l,oo). To 
account for complete absorption at the boundaries we 
let P(0,E,t) = P(L,E,t) = 0. Moreover, P(x,l,t) = 
f(x,t) is determined as a solution of 



/ = Af"-7/ + 7 



(3) 



where we approximated s» 7. This condition may 
appear somehow unusual and is justified as follows: the 
number of particles with unit energy increases at rate 7 
and are free to diffuse without increasing their energy. 
Their number diminishes by a term —7/ because they 
gain energy by amplification. Note that Eq. ((3]) can be 
derived exactly from the underlying master equation for 
the discrete variables [25j ]. 

The stationary solution of Eq. @ can be found by sep- 
aration of variables P(x, E) = Q(x)W(E): yielding Q ~ 
sin(fcc) with k = mir/L (m integer) being the separation 
constant (the wavenumber) that somehow couples diffu- 
sion and gain. The equations for W may be integrated 
exactly, but for our purposes it suffices to solve an ap- 
proximated form where energy-diffusion term in Eq. ([2]) 

is neglected, W « exp ^—Dk 2 J E T^ 1 (E')dE'^j /T(E). 

Thus for the linear gain W has a power-law tail while 
for the saturating case 



W{E) cx 



exp(-Dk 2 E/~/E s 



(4) 



(we have also assumed E s ^> E 3> 1). The general 
solution is a sum over the allowed ks. As a further 
approximation, we consider only the first Fourier mode 
k = tt/L, which is mathematically justified noting that 
Q is fixed by the stationary solution of Eq. ^ with 
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is (0 = VtTd) 

sinh(/3(ir, — L)) — sinh(/3a;) 



f + 



sinh(/3L) 



(5) 



In the critical region (to be defined below), (3 ~ tt/L so 
that / (and thus Q) is indeed very close to the shape 
of the first Fourier mode itself. The result is that the 
distribution of energies displays a parameter-dependent 
power-law 



P ( ^)ocsin(f) ^-f/^ 



(6) 



where 7 C = D{tt/L) 2 and a = 7 c /7 which is exponen- 
tially cut-off at E ~ E s . The origin of the power-law 
([6|) is traced back to very long paths which, although 
exponentially rare, acquire an exponentially large energy 
while diffusing throughout the sample 

(HE!]. Note that 

7 = 7 C correspond to the case of a Cauchy-like tail a = 1, 
P{x, E) cx E~ 2 . yielding a diverging average value of the 
energy for E s — > 00. It is thus natural to identify this 
as the "laser threshold" for the model: it will be shown 
below that this coincide with the usual criterion of gain 
overcoming the losses. 

Sofar we have described the system in terms of a sin- 
gle walker properties. Suppose we are dealing from a 
population of M walkers, and denote by Xi and Ei their 
position and energies, each obeying the equation of mo- 
tion H]). We thus describe the model in terms of the 
energy density field 



A/(t) 

4>(x,t) = EiS(x- Xi (t)) 



(7) 



The number M(t) fluctuates in time due to the fact that 
walkers are created (at a rate 7L) and absorbed at the 
boundaries (at a rate D(tt/L) 2 M), so that M ~ jL 3 /D. 
For 7 ~ 7 C ~ D I L 2 (which is the regime of interest here) 
M/L ~ I i.e. there are a few walkers per unit length and 
fluctuations must be very relevant there. To study the 
problem from this point of view, we derive the Langevin 
equation for the field. This is accomplished by standard 
Ito calculus [26j | upon writing down the equation for </>(i+ 
At) — (j)(t), expanding definition (O to second order in At 
and using Eq. 0} along with some properties of the Dirac 
S— functions [25|. If, as above, we neglect the fluctuating 
term proportional to %/T in Eq.([T]) the calculation yields 



Tt =D dx- 2+m + Tx {vq 



(8) 



The same T as in the single-walker case is used, which is 
justified since the walkers are non-interacting. The addi- 
tive, spatially uncorrelated, spontaneous emission noise 
s(x,t) is a Poissonian process whereby <f> is increased by 
e at random times, whose separation r is a distributed 
as 7exp(— 77-). The v(x,t) is a Gaussian, uncorrelated 
noise 



(v(x, t)v(x',t')) = D£S{x - x')5(t - t') 



(9) 



where the spatial scale I is introduced for dimensional 
consistency. 
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The mean-field equation obtained by neglecting fluc- 
tuations in Eq. JSJ and replacing s by its average 7 is ba- 
sically a simplified, one-dimensional version of Letokhov 
equation for random lasers Q. The steady state solu- 
tion 4>{x) is not identically zero due to the 7 term. If we 
neglect this, we straightforwardly obtain that 0, desta- 
bilizes for 7 = 7 C thus defining a threshold as assumed 
above [nj]. 

As usual for stochastic partial differential equation [29| 
Eq. ([8]) is intended as a limit of some discretization on 
a finite mesh whose spacing we denote by Ax. For defi- 
niteness and for actual numerical investigation, we choose 
the discretized equation for 0j to be (i — 0, . . . , N + 1, 
(N + l)Ax = L) 

0* = D[</. i+1 + ^_i-2^]+r(^) (10) 

with Vi being independent Gaussian variables with (vi) — 

and (vf) — D£ (we set the mesh spacing Ax = 1). The 
boundary conditions are 0o = 4>n+i = and we also im- 
pose v\ = vn — to ensure that the multiplicative term 
conserves the total energy £\ 0j. For the time derivative, 
we use a simple Euler discretization with a time step At. 

In the simulations, we choose Si to be either a Poisson 
process, (i.e. Si — 1 or zero with probability ^yAt and 

1 — 7 At, respectively) or a Gaussian with the same av- 
erage, Si — 7 + y/j^i, n being normally-distributed and 
independent random numbers with zero average and unit 
variance. Although the first choice is closer to the origi- 
nal formulation of the model, we found that in practice, 
the two processes yield almost indistinguishable results. 




FIG. 1: (Color online) Intermittent evolution of the intensity 
4>{x,t) obtained by numerical integration of Eq. (|10[l ; L — 16, 
7 = 1.027c, D = 1/2, £ = 1, r(0) = 70/(1 + 0/0*), 0* = 10 s ; 
here and in the following At — 0.01. Additive noise Sj is 
Poissonian (see text). 

The term in v of Eq. © is the leading stochastic cor- 
rection to the mean-field evolution. It is a multiplica- 



tive process and can be regarded as a kind of random 
advection whereby fluctuations are transported almost 
coherently, while keeping the total energy conserved. It 
thus acts against the diffusive term that tends to smooth 
out fluctuations. As a consensequence, the field is highly 
intermittent in time, with large-amplitude bursts emerg- 
ing from a lower-amplitude background (Fig. [T]). It is 
thus intuitively plausible that such term is responsible of 
most of the nontrivial statistics of the field. As a mat- 
ter of fact, random advection dynamics is known to yield 
strongly non-Gaussian distributions for the fields, and 
even power-law tails [3(| . Actually, equations similar to 
(|10p (with linear gain and periodic boundary conditions) 
have been studied in Refs. [3lj, 32] as a model for an active 
scalar (e.g. a temperature field) convected by a random 
velocity field. It was argued that power law tails generi- 
cally arise when both multiplicative and additive noises 
are present, which is indeed the case here. As a matter 
of fact, there are two sources of additive noise. One is of 
course the spontaneous emission term s; the other stems 
from the fact that the deterministic, steady state value 
is non vanishing, thus yielding a additive contribution 
of order d(v(f>)/dx. 

There is however a crucial difference, which is dictated 
by the physical origin of the noise: the advective term is 
a finite-size effect which decreases with the system size. 
This can be demonstrated by dimensional analysis. To 
be more general, let us consider the extension of Eq. (|SJ) 
to d dimensions, 

^ = £V 2 + 70 + V-(v0)+ 7 (11) 
at 

where v is a d— dimensional vector whose components 
v v (y — 1 , . . . , d) are Gaussian distributed and satisfy a 
relation akin to ©, 

{v v (x,t)v v/ (x',t')) = D£ d 8 v ^8 d (x-x')8{t-t') (12) 

For simplicity, we neglected again the VT term, replaced 
s by its average 7 and restricted to the simpler case 
P(0) = 70. Let us rewrite Eq. (|8|) introducing the di- 
mensionless variables x/L —> x, Dt/L 2 — > t. <p/e — > 
and we also rescale Dj/L 2 — > 7. The noise rescales ac- 
cordingly as yj L d+2 j D v — > v. Thus, Eq. (J8)) is written 
in the dimensionless form 

^ = V 2 + 70 + A(L)V-(u0) + 7 (13) 

where u is again an uncorrelated Gaussian noise with 
unit variance and we have made explicit the rescaled 
noise strength which turns out to be A cx (i/L)?. As 
a consequence, we expect that the multiplicative noise to 
yield sizeable power-law tails only for small enough sizes. 
When L — V 00, the diffusive term in Eq. (fTT|) dominates 
and power law are no longer observable. 

The numerical simulations of Eq. (|10[) are in qualitative 
agreement with the above argument. We monitored the 
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distribution of the field at the center of the mesh 4>n/2 a s 
well as the the outcoming flux J = D[<f>x(t) + <p^(t)]/2 as 
function of time. This quantity is of experimental inter- 
est, being related to the emission spectra. The tail of the 
distributions of both observables display, upon increas- 
ing L, a smooth crossover from a power-law tail with an 
exponent a very close to the one predicted by Eq. ([6]) to 
a faster decay (see FigE^i). The same occurs for fixed L 
upon approaching the threshold j c (FigEb). 

The existence of fat tails is inimately related to the 
the possibility for a spontaneous fluctuation to grow well 
beyond the average. The indicator to quantify this is the 
generalized Lyapunov exponent 33(. For a perturbation 
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8<j){x, t) of the field, which evolves according to the lin- 
earized version of Eq.®, let R(t) = \\8cf>(t + T)\\/\\5(j>{t)\\ 
be the response function after a time r to a disturbance 
at time t. The generalized finite-time Lyapunov expo- 
nent A(q) is then defined by R q (r) ~ exp(A(g)r) where 
the overline denote an average over time. If A(q) > for 
large enough q then the system has a finite probability 
that a small perturbation results in a large fluctuation. 
Moreover, the deviations of A(q) from a linear behaviour 
in q are a measure of intermittency [33| . Figj3] shows 
that A becomes positive for larger and larger values of 
q upon increasing L, signalling that events far from the 
average becomes increasingly rare, in agreement with the 
dimensional arguments discussed above. 




FIG. 2: (Color online) Probability distribution functions of 
the flux J for (a) different L, 7 = 7 C and (b) L — 16, different 
7 values. The straight dashed line corresponds to the power 
law a — 1 expected from Eq. (0 . Flux data have been binned 
over consecutive time windows of duration Tw = 10, over a 
run of about 10 6 time units. Additive noise is Gaussian (see 
text); other parameters as in previous figure. 



To summarise, we have presented a simple model for 
diffusive random media with gain which yields power law 
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FIG. 3: (Color online) The generalized finite-time Lyapunov 
exponent A(q) for r = 200 and different system lengths L. 



distribution of the intensities. Our main result is the 
Langevin equation for the energy density field, Eq.®, 
that establishes a novel connection between the physics 
of scattering media with gain (e.g. random lasers) and 
the theory of noncquilibrium phenomena in spatially ex- 
tended systems. The random advection, multiplicative 
term in Eq.® competes with diffusive and gain terms 
and is responsible for unconventional fluctuations. Its 
relevance is gauged by the effective noise strength X(L). 
This means that Levy-like fluctuations are more likely to 
be observed in lower dimensions (for instance in d = 2) 
and when the mean free path of photons is not too short 
with respect to the sample size. Those simple criteria 
should be of general guidance to interpret and compare 
data for different type of scattering media and in different 
experimental condition. 

I thank S. Cavalieri, F. Ginelli and R. Livi for useful 
comments. 
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